Experimental enforcement of sex-limited adaptation

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, cache = FALSE)

Setup for analysis

Load packages

library(tidyverse) # tidy style coding
library(brms) # Bayesian models
library(loo) # to use information criteria in brms models
library(tidybayes) # Bayesian aesthetics
library(MetBrewer) # colours
library(kableExtra) # tables
library(patchwork) # putting plots together
library(DT) # for search- and saveable tables
library(pander) # for simpler tables
library(ggdag) # to visualise casual relationships

Load in the data

\(~\)

fitness_data <- read_csv("data/SLC_fitness_data.csv") %>% 
  mutate(Fitness_vial_ID = as.factor(Fitness_vial_ID),
         Block = as.factor(Block),
         Population = as.factor(Population),
         Treatment = as.factor(Treatment),
         GFP = as.factor(GFP),
         Sex = as.factor(Sex),
         Rearing_vial = as.factor(Rearing_vial),
         Sex = as.factor(Sex),
         Total_red_offspring = Red_female_offspring + Red_male_offspring,
         Total_bw_offspring = bw_female_offspring + bw_male_offspring,
         Total_offspring = Total_red_offspring + Total_bw_offspring) %>% 
  rename(Evolution_treatment = Treatment)

# Create a function to build HTML searchable tables

my_data_table <- function(df){
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons =
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'fitness_data')),
      pageLength = 692
    )
  )
}


my_data_table(fitness_data %>% select(-Comment))

Column explanations

Fitness_vial_ID: a unique identifier for each replicate of the fitness assays.

Block: the experiment was run in three distinct blocks, using flies from separate generations.

Population: we measured the fitness of flies from 12 independent populations that has undergone experimental evolution.

Evolution_treatment: the populations had been exposed to one of three evolutionary conditions for 20 generations: a female-limited response to selection, a male-limited response and a control condition where the evolutionary response was unconstrained.

GFP: the GFP marker carried by the population. UBI indicates the presence of a transgene that encodes ubiquitous expression of GFP, while 3xP indicates the presence of a different transgene that encodes the expression of GFP in the ocelli.

Sex: the sex of the individuals that we were measuring the fitness of.

Rearing_vial: the vial the flies developed in. This variable is included to capture variation explained by the rearing environment e.g. small differences in food moisture content or quantity. Note that females and males can have the same rearing vial.

Red_female_offspring: the number of adult female offspring sired/produced by flies sourced from one of the 12 populations.

Red_male_offspring: the number of adult male offspring sired/produced by flies sourced from one of the 12 populations.

bw_female_offspring: the number of adult female offspring sired/produced by the competitor flies in our fitness assay. bw is a recessive allele that encodes brown eye colour.

bw_male_offspring: the number of adult male offspring sired/produced by the competitor flies in our fitness assay. bw is a recessive allele that encodes brown eye colour.

Total_red_offspring: the total number (sexes pooled) of adult offspring sired/produced by flies sourced from one of the 12 populations.

Total_bw_offspring: the total number (sexes pooled) of adult offspring sired/produced by competitor flies.

Total_offspring: the total number (sexes and eye colours pooled) of adult offspring counted in each vial.

\(~\)

Modelling approach

\(~\)

Female and male fitness are fundamentally different concepts / traits. There are also several differences between our female and male fitness assays. The major difference is that the male assay contains half the number of females in any given vial than does the female assay. The logic behind this design choice is that sexually selected processes are a more important determinant of male fitness than they are female fitness, so any fitness differences may only be observed when competition for fertilisations is high.

For these reasons, we choose to split the data up and model them separately.

female_fitness <- 
  fitness_data %>%
  filter(Sex == "Female")

male_fitness <- 
  fitness_data %>%
  filter(Sex == "Male")

Fixed effects

Block: fitness might differ between the three distinct blocks we split our experiment up into. Blocks differed temporally, used flies from different generations and different batches of food. It is also possible that there were minor fluctuations in the lighting and temperature environment experienced during development between blocks. Each of these variables may introduce variation into our fitness measurements, that can be accounted for by including the Block variable in our model.

GFP: it is possible that fitness may be affected by the GFP transgene carried by each population. For example, one could imagine that any unintended fitness effects of a transgene might be of greater magnitude if it is expressed in a larger proportion of tissues, as is the case for the UBI transgene versus the 3xP transgene. Note that each GFP type is carried by an equal number of populations from each of the three evolutionary treatments.

Evolution_treatment: this is the evolution condition that populations were subject to for 20 generations. There are three levels: populations carrying female-adapted chromosomes, populations containing male-adapted chromosomes and populations carrying control chromosomes that have experienced an unlimited response to selection. This variable is central to our hypothesis.

Random effects

Rearing_vial: the vial individual flies developed within may introduce further variation into our response variable. Like Block this variable controls for micro-environmental variation.

Population: our design contained 12 independent populations of chromosomes that originated from a single outbred laboratory fly population. The populations were split and chromosomes from each were subjected to one of the three evolution treatments for 20 generations. 4 populations experienced a female-limited response to selection, 4 a male-limited response and 4 an unlimited response.

\(~\)

Female fitness

\(~\)

To estimate the fitness of females carrying each of the three chromosome types, we placed three experimental females into a yeasted vial with three female competitors that carried the bw mutation. We then introduced six males that also carried the bw mutation. We allowed them to mate and oviposit for three days, then removed all adults from the vial. 12 days later we counted all of the adult progeny in the vial and scored them for eye-colour. Progeny with red eyes were produced by the experimental females ( bw is recessive) and progeny with brown eyes were produced by the competitor females. We calculated fitness as the proportion of red eyed offspring in the vial.

\(~\)

female_fitness.01 <-
  brm(Total_red_offspring | trials(Total_offspring) ~ 1 + Evolution_treatment + GFP + Block + (1|Population) + (1|Rearing_vial),
      data = female_fitness, family = binomial(), 
      prior = c(prior(normal(0, 1.5), class = Intercept),
                prior(normal(0, 2), class = b),
                prior(exponential(1), class = sd)),
      iter = 10000, warmup = 8000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.999, max_treedepth = 15),
      seed = 2, file = "Fits/female_fitness.01")

#add_criterion(female_fitness.02, criterion = "waic")

print(female_fitness.01)

##  Family: binomial 
##   Links: mu = logit 
## Formula: Total_red_offspring | trials(Total_offspring) ~ 1 + Evolution_treatment + GFP + Block + (1 | Population) + (1 | Rearing_vial) 
##    Data: female_fitness (Number of observations: 327) 
##   Draws: 4 chains, each with iter = 10000; warmup = 8000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Population (Number of levels: 12) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.06      0.05     0.00     0.19 1.01      649     1437
## 
## ~Rearing_vial (Number of levels: 70) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.27      0.03     0.22     0.33 1.00     1772     3026
## 
## Population-Level Effects: 
##                                   Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                             1.09      0.10     0.91     1.28 1.00
## Evolution_treatmentFemale_limited     0.19      0.10    -0.01     0.39 1.00
## Evolution_treatmentMale_limited       0.16      0.10    -0.03     0.36 1.00
## GFPUBI                               -0.32      0.08    -0.48    -0.16 1.00
## Block2                               -0.58      0.08    -0.75    -0.43 1.00
## Block3                               -0.61      0.08    -0.77    -0.45 1.00
##                                   Bulk_ESS Tail_ESS
## Intercept                             2167     3294
## Evolution_treatmentFemale_limited     2313     2741
## Evolution_treatmentMale_limited       2397     3541
## GFPUBI                                2188     3499
## Block2                                1879     3162
## Block3                                1882     3143
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Conduct a posterior predictive check to confirm our model is doing what we want it to.

pp_check(female_fitness.01, type = "hist", ndraws = 11, binwidth = 10) +
  theme_minimal() +
  theme(panel.background = element_blank())

The posterior predictive distribution matches our raw data quite well, indicating the model is functioning as we wanted.

\(~\)

Get predictions from the model and present them in a Table

draws <- as_draws_df(female_fitness.01)

draws_female <-
  draws  %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Evolution_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Evolution_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>% 
  select(Control, `Female-limited`, `Male-limited`) %>% 
    pivot_longer(cols = c(Control, `Female-limited`, `Male-limited`),
                 names_to = "Evolution_treatment") %>% 
  rename(prop_focal_offspring = value)

draws_female %>% 
  group_by(Evolution_treatment) %>% 
  summarise(`Estimated prop. of offspring produced` = median(prop_focal_offspring),
            `2.5%` = quantile(prop_focal_offspring, probs = 0.025),
            `97.5%` = quantile(prop_focal_offspring, probs = 0.975)) %>% 
  rename(`Evolution treatment` = Evolution_treatment) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Evolution treatment	Estimated prop. of offspring produced	2.5%	97.5%
Control	0.749	0.712	0.783
Female-limited	0.784	0.749	0.815
Male-limited	0.778	0.746	0.808

# now find the differences between the control and the sex-limited treatments

# the inv_logit_scaled() function converts the posterior draws onto the response scale 

  draws %>% 
  mutate(p_control =  inv_logit_scaled(b_Intercept),
         p_female = inv_logit_scaled(b_Evolution_treatmentFemale_limited + b_Intercept),
         p_male = inv_logit_scaled(b_Evolution_treatmentMale_limited + b_Intercept),
         `Female-limited` = p_female / p_control,
         `Male-limited` = p_male / p_control) %>% 
  gather(key = `difference comparison`, value = `% difference`) %>% 
  filter(`difference comparison` == c("Female-limited", "Male-limited")) %>% 
  group_by(`difference comparison`)  %>% 
  summarise(`Mean proportion of offspring produced relative to control`  = mean(`% difference`),
            `2.5%` = quantile(`% difference`, probs = 0.025),
            `97.5%` = quantile(`% difference`, probs = 0.975)) %>% 
  rename(`Evolution treatment` = `difference comparison`) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Evolution treatment	Mean proportion of offspring produced relative to control	2.5%	97.5%
Female-limited	1.046	0.997	1.097
Male-limited	1.039	0.993	1.091

\(~\)

Male fitness

\(~\)

To estimate the fitness of males carrying each of the three chromosome types, we conducted an experiment very similar to the female fitness assay. However, because male fitness has stronger covariance with fertilisation events than does female fitness, we conducted the male fitness assay with a 1:2 sex ratio (female:male) rather than the 1:1 ratio used in the female assay. This increases the strength of sexual selection and is a more appopriate way to expose differences in fitness between groups of males. As with the females, we calculated fitness as the proportion of red eyed offspring in the vial.

male_fitness.01 <-
  brm(Total_red_offspring | trials(Total_offspring) ~ 1 + Evolution_treatment + GFP + Block + (1|Population) + (1|Rearing_vial),
      data = male_fitness, family = binomial(), 
      prior = c(prior(normal(0, 1.5), class = Intercept),
                prior(normal(0, 2), class = b),
                prior(exponential(1), class = sd)),
      iter = 10000, warmup = 8000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.999, max_treedepth = 15),
      seed = 2, file = "Fits/male_fitness.01")

#add_criterion(male_fitness.02, criterion = "waic")

print(male_fitness.01)

##  Family: binomial 
##   Links: mu = logit 
## Formula: Total_red_offspring | trials(Total_offspring) ~ 1 + Evolution_treatment + GFP + Block + (1 | Population) + (1 | Rearing_vial) 
##    Data: male_fitness (Number of observations: 360) 
##   Draws: 4 chains, each with iter = 10000; warmup = 8000; thin = 1;
##          total post-warmup draws = 8000
## 
## Group-Level Effects: 
## ~Population (Number of levels: 12) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.29      0.15     0.05     0.63 1.00     1284     1796
## 
## ~Rearing_vial (Number of levels: 72) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.62      0.07     0.51     0.76 1.00     2655     4566
## 
## Population-Level Effects: 
##                                   Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept                             0.66      0.26     0.16     1.16 1.00
## Evolution_treatmentFemale_limited     0.13      0.29    -0.43     0.69 1.00
## Evolution_treatmentMale_limited       0.27      0.29    -0.33     0.82 1.00
## GFPUBI                               -0.07      0.23    -0.54     0.40 1.00
## Block2                                0.92      0.18     0.56     1.28 1.00
## Block3                                0.05      0.18    -0.30     0.39 1.00
##                                   Bulk_ESS Tail_ESS
## Intercept                             5450     5236
## Evolution_treatmentFemale_limited     5534     4821
## Evolution_treatmentMale_limited       5868     5011
## GFPUBI                                6426     5062
## Block2                                4454     5090
## Block3                                4342     4641
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Conduct the posterior predictive check…

pp_check(male_fitness.01, type = "hist", ndraws = 11, binwidth = 10) +
  theme_minimal() +
  theme(panel.background = element_blank())

Get predictions from the model and present them in a Table

draws_m <- as_draws_df(male_fitness.01)

draws_male <-
  draws_m  %>% 
  mutate(`Female-limited` = inv_logit_scaled(b_Intercept + b_Evolution_treatmentFemale_limited),
         `Male-limited` = inv_logit_scaled(b_Intercept + b_Evolution_treatmentMale_limited),
         Control = inv_logit_scaled(b_Intercept)) %>% 
  select(Control, `Female-limited`, `Male-limited`) %>% 
    pivot_longer(cols = c(Control, `Female-limited`, `Male-limited`),
                 names_to = "Evolution_treatment") %>% 
  rename(prop_focal_offspring = value)

draws_male %>% 
  group_by(Evolution_treatment) %>% 
  summarise(`Estimated prop. of offspring sired` = median(prop_focal_offspring),
            `2.5%` = quantile(prop_focal_offspring, probs = 0.025),
            `97.5%` = quantile(prop_focal_offspring, probs = 0.975)) %>% 
  rename(`Evolution treatment` = Evolution_treatment) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Evolution treatment	Estimated prop. of offspring sired	2.5%	97.5%
Control	0.659	0.54	0.761
Female-limited	0.686	0.576	0.782
Male-limited	0.716	0.6	0.807

# now find the differences between the control and the sex-limited Evolution_treatments

# the inv_logit_scaled() function converts the posterior draws onto the response scale 

draws_m %>% 
  mutate(p_control =  inv_logit_scaled(b_Intercept),
         p_female = inv_logit_scaled(b_Evolution_treatmentFemale_limited + b_Intercept),
         p_male = inv_logit_scaled(b_Evolution_treatmentMale_limited + b_Intercept),
         `Female-limited` = p_female / p_control,
         `Male-limited` = p_male / p_control) %>% 
  gather(key = `difference comparison`, value = `% difference`) %>% 
  filter(`difference comparison` == c("Female-limited", "Male-limited")) %>% 
  group_by(`difference comparison`)  %>% 
  summarise(`Mean proportion of offspring sired relative to control`  = mean(`% difference`),
            `2.5%` = quantile(`% difference`, probs = 0.025),
            `97.5%` = quantile(`% difference`, probs = 0.975)) %>% 
  rename(`Evolution treatment` = `difference comparison`) %>% 
  pander(split.cell = 40, split.table = Inf, round = 3)

Evolution treatment	Mean proportion of offspring sired relative to control	2.5%	97.5%
Female-limited	1.047	0.872	1.255
Male-limited	1.093	0.906	1.304

Building Figure 1

# female plots

f_1 <-
  draws_female %>% 
  ggplot(aes(Evolution_treatment, prop_focal_offspring)) + 
  stat_eye(aes(fill = Evolution_treatment), .width = c(0.66, 0.95), alpha = 1) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
  scale_y_continuous("Female fitness\n(proportion of offspring produced)",
                     breaks = c(0.4, 0.6, 0.8)) +
  coord_cartesian(ylim = c(0.4, 1)) +
  xlab(NULL) +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

f_2 <-
  draws %>% 
  mutate(`Female-limited` = b_Evolution_treatmentFemale_limited,
         `Male-limited` = b_Evolution_treatmentMale_limited) %>% 
  gather(key = parameter, value = logodds) %>% 
  filter(parameter == c("Female-limited", "Male-limited")) %>%
  as_tibble() %>% 
  
  ggplot(aes(parameter, logodds)) + 
  stat_halfeye(aes(fill = parameter), .width = c(0.66, 0.95)) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = c("Female-limited" = "#ffd06f", "Male-limited" = "#72bcd5")) + 
  coord_flip(ylim = c(-1, 1)) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_y_continuous(breaks = c(-1, 0, 1)) +
  xlab(NULL) +
  ylab(NULL) +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

# male plots

f_3 <-
  draws_male %>% 
  ggplot(aes(Evolution_treatment, prop_focal_offspring)) + 
  stat_eye(aes(fill = Evolution_treatment), .width = c(0.66, 0.95), alpha = 1) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = met.brewer("Hiroshige", 3)) +
  scale_y_continuous("Male fitness\n(proportion of offspring sired)",
                     breaks = c(0.4, 0.6, 0.8)) +
  coord_cartesian(ylim = c(0.4, 1)) +
  xlab("Evolutionary history") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())

f_4 <-
  draws_m %>% 
  mutate(`Female-limited` = b_Evolution_treatmentFemale_limited,
         `Male-limited` = b_Evolution_treatmentMale_limited) %>% 
  gather(key = parameter, value = logodds) %>% 
  filter(parameter == c("Female-limited", "Male-limited")) %>%
  as_tibble() %>% 
  #mutate(parameter =factor(parameter, levels=c("SD-Mad", "SD-72", "SD-5"))) %>%
  
  ggplot(aes(parameter, logodds)) + 
  stat_halfeye(aes(fill = parameter), .width = c(0.66, 0.95)) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals
  scale_fill_manual(values = c("Female-limited" = "#ffd06f", "Male-limited" = "#72bcd5")) + 
  coord_flip(ylim = c(-1, 1)) +
  geom_hline(yintercept = 0, linetype = 2) +
  scale_y_continuous(breaks = c(-1, 0, 1)) +
  xlab(NULL) +
  ylab("Log-odds mean difference from control") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank())


(f_1 + f_2) /(f_3 + f_4)

Figure 1: a shows the estimated mean female fitness for flies carrying autosomes that had previously experienced unconstrained inheritance (control), female-limited inheritance or male-limited inheritance. b shows the log-odds mean difference in means between the sex-limited treatments and the control evolutionary treatment. c and d depict the same things as a and b, except for estimated mean male fitness.